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Abstract. In this paper we present and discuss results of Monte Carlo numerical simulations of the 
two-dimensional Ising ferromagnet in contact with a heat bath that intrinsically has a thermal gradient. 
The extremes of the magnet are at temperatures Ti < T c < T2, where T c is the Onsager critical 
temperature. In this way one can observe a phase transition between an ordered phase (T < T c ) and a 
disordered one (T > T c ) by means of a single simulation. 

By starting the simulations with fully disordered initial configurations with magnetization m = 
corresponding to T = 00, which are then suddenly annealed to a preset thermal gradient, we study the 
short-time critical dynamic behavior of the system. Also, by setting a small initial magnetization m — mo, 
we study the critical initial increase of the order parameter. Furthermore, by starting the simulations 
from fully ordered configurations, which correspond to the ground state at T — and are subsequently 
quenched to a preset gradient, we study the critical relaxation dynamics of the system. Additionally, we 
perform stationary measurements (t —¥ 00) that are discussed in terms of the standard finite-size scaling 
theory. 

We conclude that our numerical simulation results of the Ising magnet in a thermal gradient, which are 
rationalized in terms of both dynamic and standard scaling arguments, are fully consistent with well 
established results obtained under equilibrium conditions. 
PACS Numbers : 64.60.De; 64.60. an, 05.10.-a 

Keywords : Phase transitions, thermal gradient, Ising ferromagnet, numerical simulations. 



PACS. XX.XX.XX No PACS code given 



2 Muglia and Albano: The Ising 

1 Introduction 

Let us assume that a physical system undergoes a smooth 
transition between two distinctly characteristic phases when 
a certain control parameter (e.g., T = temperature, P = 
pressure, fi — chemical potential, etc.) is finely tuned 
around a critical value. This physical situation is often 
addressed by means of analitical and numerical methods 
where the control parameter asuumes a fixed value. How- 
ever, here we considered an alternative approach by as- 
suming that along a given spatial direction, the system 
undergoes a well established gradient of the control pa- 
rameter, such that the critical point lies between the ex- 
treme values, e.g., T\ < T c < T2, where one has a ther- 
mal gradient between temperatures T\ and T2, and the 
critical temperature T c is within that range [U^EHEj • 
Of course, a standard reversible system, usually studied 
under equilibrium conditions, will be out of equilibrium 
under the gradient constraint applied to the control pa- 
rameter. However, we will show below that this situation 
is no longer a shortcoming for the application of methods 
and theories already developed for the study of equilib- 
rium systems. In order to fix ideas, in this paper we will 
study the dynamical and stationary behavior of the Ising 
model for a two-dimensional magnet 6 in a thermal gra- 
dient. The proposed study not only possess interesting 
theoretical challenges, but it may also be useful in con- 
nection to recent experimental work aimed to characterize 
films in general, and magnetic films in particular, which 
are obtained under thermal gradient conditions imposed 
on the substrate. In fact, since the temperature is a key 
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parameter for the relevant properties of thin films, several 
experiments have focused on the influence of a tempera- 
ture gradient on film growth, e.g., Tanaka et al. reported 
studies of magnetic Tb-Te films obtained under thermal 
gradient conditions [7] ■ Also, Schwickert et al. |5] have in- 
troduced the temperature wedge method where a temper- 
ature gradient of several hundred Kelvin was established 
across a substrate during the co-deposition of Te and Pt. 
Other experiments performed by Xiong et al.[5] involve 
nanostructures obtained by using a gradient temperature. 

So, within the broad context discussed above, we now 
focus our attention on the aim of this paper, which is to 
perform an extensive numerical study of the two-dimensional 
Ising ferromagnet in a thermal gradient. In previous nu- 
merical studies performed by various authors [U^EIHE] , 
the two extremes of the magnet are considered in con- 
tact with different thermal baths, e.g. the left-hand and 
the right-hand side extremes are at fixed temperatures T\ 
and Ti , respectively. In order to study this particular out- 
of-equilibrium situation, varoius algorithms can be used, 
e.g. the Creutz algorithm [TlllOj. the Kadanoff-Swift move 
[TT] the Q2R rule [12], the KQ rule [3], a recently intro- 
duced microcanonical dynamics [5], etc. In this way, the 
system naturally evolves into a stationary state and a ther- 
mal gradient is established between its extremes. Also, the 
surfaces in the direction perpendicular to the one where 
the gradient is observed are surrounded by an adiabatic 
wall. In particular, Neek-Amal et al. [2j have reported the 
ocurrence of an almost linear temperature gradient (c.f. 
figure 1 of reference [5]), for the case of the d — 2 Ising 
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model, by using the Creutz algorithm. On the other hand, 
the mean-field treatment given by the well known Fourier 
equation [13] . namely 

d 2 T(x,t) _ K dT(x,t) = Q 
dx 2 dt 

T(0,i) = Ti, T(L,t)=T 2 , 

where T(x, t) is the temperature along the solid and K 
is the (constant) thermal conductivity, predicts a linear 
temperature gradient between both baths. In fact, we as- 
sumed a solid of length L where the gradient is applied, 
while the remaining directions are irrelevant. Then, in the 
t — > oo limit, the time-dependent contributions to the so- 
lution become negligible, and one gets the linear gradient, 
namely, 

T{x) =T 1+ {T2 ~ Tl) x. (1) 

It is worth mentioning that the above simulation stud- 
ies [TJ[2l[3l|5] of the Ising system are precisely focused on 
the calculation of the thermal conductivity, among other 
relevant observables, such as the energy profiles of the sys- 
tems transversal sections. Also, in all cases the discussed 
results were obtained under stationary conditions. How- 
ever, simulation results that are in contrast to the mean- 
field theory, clearly show that the thermal conductivity is 
not a constant in these systems since it depends on the 
local temperature (or energy) QQE212Jt2L bence the tem- 
prature may not necessarily grow linearly along the heat 
propagation direction. In fact, in the presence to ther- 
mal fluctuations it is expected that transport properties 
may exhibit non-trivial spatial patterns due to the non- 
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equilibrium nature of the system under study [5] , and un- 
der these conditions the thermal conductivity depends on 
the temperature (apart from other relevant parameters of 
the model). In particular, for the Ising system in a ther- 
mal gradient, Agliari et al. [5] reports that K(T) exhibits 
a peak slightly above criticality, whose position is inde- 
pendent of both the lattice size and the temperature dif- 
ference between the extrems of the sample. Therefore, the 
system that we study in the present paper is not simply 
an Ising model where different temperatures arc imposed 
at the bundaries, but it describes a system where the tem- 
perature gradient is present in the heat bath itself and the 
transport properies of the system are determined also by 
the bath and not only by the Ising dynamics. While the 
linear gradient used in our calculations is motivated by the 
fact that we considered a sample in contact with a thermal 
bath that intrinsically has a gradient, that assumption can 
also be supported by a physical argument. In fact, to de- 
fine a thermal conductivity there must exist mechanisms 
such that the degrees of freedom (e.g. the phonons) dom- 
inate the thermal properties and spins locally thermalize 
to the temperature of the lattice. In the absence of such 
mechanisms the phonons at one end of the crystal will not 
be in thermal equilibrium at a temperature T 2 , and those 
at the other end in equilibrium at T\ |13) . Then, if we as- 
sume that the lattice has a constant conductivity, at least 
in the temperature range where the magnetic transition 
takes place, we expect a linear growth of the temperature 
when a temperature difference is apllied to the extrems. 
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In the present paper we not only undertake a com- fective dimension d* = 1, which reflects the fact that the 
pletely different approach than those used in previous stud- susceptibility is measured as the fluctuations of the order 
ies, as already anticipated, but we are also interested in parameter (magnetization) per unit of length along the 
both the dynamic and the stationary behavior of a broad direction perpendicular to the gradient, 
spectra of physical observables. In fact, we consider that The paper is organized such that in Section [5] we pro- 
the Ising ferromagnet is already in contact with a ther- vide the description of the model and the simulation method, 
mal bath that intrinsically has a thermal gradient along in order to allow for a smooth introduction of the theoret- 
the x— direction, as described by equation ((!]). In this way ical background in Section [31 In section U] we present and 
each column of coordinates Xi — i, with i — 1,2, ....,L discuss the results. Finally, our conclusions are stated in 
of our two-dimensional (lattice) system, is in equilibrium Section [5j 
with the thermal bath at temperature T(xi). Under this 
condition we apply the standard Metropolis dynamics by 
considering the proper temperature for each column. So, 
the whole system reaches a stationary, non equilibrium 
regime, with a net flux of energy from the hotter to the 
colder extremes at temperatures T2 and T±, respectively. 



2 BRIEF DESCRIPTION OF THE ISING 
MAGNET IN A THERMAL GRADIENT 
AND THE SIMULATION METHOD 



We performed simulations of the Ising model in the two- 
As already anticipated, we are interested in the study dimensional square lattice with a rectangular geometry, 
of the dynamic behavior of the system, i.e., by addressing with a Hamiltonian given by 
both the relaxation dynamics and the short-time dvnamics|14|. 

H = -J SijSi'j', (2) 

as well as the stationary behavior that is established in 

the long-time regime. By using the dynamic scaling the- where the spin variables can assume two values, i.e. Sij = 
ory [14] and the finite-size scaling theory [18], both de- ±1, J > is the coupling constant that is set positive 
veloped for the analysis of critical behavior of samples in in order to account for ferromagnetic interactions, and 
homogeneous thermal baths, we are able to rationalize all summation runs over nearest-neighbor sites as indicated 
measurements performed in our thermal gradient system, by the symbol The Ising system is an archety- 

In this way we determined not only the critical tempera- pal model for the study of second-order equilibrium phase 
ture but also relevant critical exponents such as those of transitions. In fact, it is well known that the local interac- 
the order parameter (f3), the susceptibility (7), the corre- tions of the Hamiltonian given by equation $2$ lead to a 
lation length (v), etc. We also show that the hyperscaling macroscopically observable transition between an ordered 
relationship given by d* — 2/3 /v — j/v holds, with an ef- ferromagnetic phase and a disordered paramagnetic one, 
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which in d — 2 dimensions takes place at the Onsager quently, the standard Monte Carlo procedure is imple- 
critical temperature T c — - — —7^ — 7=r ll5l . where fee is the mented. 

ks ln(l+V2) 

Boltzmann constant (hereafter we report all temperature The upper panel of figure 1 shows a typical snapshot 
values in units of J/fcs). configuration obtained under stationary conditions by us- 



In this work we placed the Ising ferromagnet of size 
L x M in contact with a thermal bath that intrinsically 
has a thermal gradient along the L— direction, as given by 
equation (fTJ) . Therefore, each column of length M located 
at x = i, with 1 < i < L, is in thermal equilibrium with 
the bath temperature T(i) = T 2 + (T 2 - T^i/L, where T x 
and Ti are the temperatures on the left- and right-hand 
sides of the magnet, respectively. Since we are interested 
in the study of the critical behavior, the temperatures are 
selected such that T\ < T c < T 2 , so that a phase transi- 
tion between the ordered phase on the left-hand side and 
the disordered one at the right-hand side can be observed 



ing the above-described procedure, while the lower panel 
shows a sketch of the linear gradient applied to the sample. 



in a single simulation run. As imposed by the simulation m[ n ^ (T(i),x,t) = ( — s i,j(t) j 



Due to the applied thermal gradient, all relevant phys- 
ical observables depend on the coordinate running in the 
direction of the gradient itself, so one essentially has to 
measure profiles of the observables, such as the nth mo- 
ment of the order parameter (magnetization in this case) 
is given by 

M 

(3) 



0- 

where the dependence on the gradient axis (1 < i < L) is 
explicitly indicated. Here, we have also included the time 
dependence, where a Monte Carlo time unit involves the 
random selection of L x M spins, as usual. 

Furthermore, other measured observables are the sus- 
ceptibility profiles, evaluated as the magnetization fluctu- 
ations, i.e., 

uated by using the Hamiltonian given by equation ©. Xi(T(i),t) = — [m[ 2 \T (i) , t) - (m[ Vj ' (T (i) , t)) 2 ]; (4) 



geometry used we take periodic boundary conditions in 
the direction perpendicular to the gradient, while open 
boundary conditions are adopted for the extremes of the 
sample in the direction parallel to the gradient. 

In order to perform the Monte Carlo simulation, a 
spin Sij selected at random and the change of the energy 
(AH) involved in a flipping attempt — > are eval- 



Now, due to the applied thermal gradient, the flipping 



and the second- and the fourth-order cumulants 



probability given by the Metropolis dynamics W(sij — > i 2 \ 

u 2 (T(i),t)= m ;J T ^ -1, (5) 

-Sij) = exp(—AH/kBT(i)) depends on the ith column (m\ (T(i),t)) 2 
(1 < i < L) where the selected spin is located. Subse- and 
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Fig. 1. (Color online) The upper panel shows a typical snap- 
shot configuration of a two-dimensional Ising ferromagnet in 
contact with a thermal bath that intrinsically has a linear tem- 
perature gradient that grows from T\ — 0.2 (left-hand side) to 
T2 = 4.338 (right-hand side). Up and down spins are shown in 
blue and red, respectively. The lower panel is a sketch of the 
applied gradient, showing the critical region around T c ~ 2.269 
by means of a dotted line. The critical interphase between the 
ordered phase (left-hand side) and the disordered one (right- 
hand side) can be observed in the snapshot. More details in 
the text. 



U 4 (T(i),t) = 1- 



3(mf(T{i),t)r 



(6) 



respectively. 

In order to understand the behavior of spin-spin cor- 
relations in a system in the presence of a temperature 



gradient, we also measured the spin-spin spatial correla- 
tion function in the direction perpendicular to the thermal 
gradient, defined for each column i at temperature Tj as 

ftf-r 



9(Ti,r) 



1 



(M-r) 



3=1 



where r is the distance between spins. Furthermore, we 
also measured the correlation function in the direction 
parallel to the thermal gradient, by fixing the origin just 
at T c . 

For the sake of completeness it is worth mentioning 
that for dynamic simulations the time dependence of the 
already defined observables is obtained by averaging over 
a (large) number of different realizations (JVr) performed 
with different sets of random numbers and physically equiv- 
alent initial conditions. Also, measurements corresponding 
to the stationary regime are performed by taking aver- 
ages after a large number (N c ) of different configurations, 
which are obtained after disregarding a (large) number 
(Nd) of Monte Carlo time steps, in order to allow for the 
system stabilization. 

3 THEORETICAL BACKGROUND 

By considering a fully disordered sample but with a small 
initial magnetization (mo), the general scaling behavior of 
the n-th moment of the magnetization that it is expected 
to hold for temperatures close to the critical region, such 
that e = (T( 'jr Tc) ~ 0, is given byQI], 

m^{e,t,L) =b-¥m(k)(b- z t,bie,b~ 1 L,b Xo m ), (8) 

where /3 and v are the order parameter and the correla- 
tion length (static) critical exponents, z is the dynamic 
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exponent, and 6 is a suitable scaling parameter [14]. Ac- while the logarithmic derivative of equation (lllj) evaluated 
cordingly to the short-time dynamic scaling theory, x$ is just at criticality gives [Tl] 



a new exponent that governs the early-time scaling be- 



d\n(m(t)) 



de 

havior of the moments of the initial magnetization. Now, 



"t 1 /"*. (13) 

e=0 



by setting b = t 1 ^, just at the critical point (e = 0[TB]), Finally, note that the scaling behavior of the cumulants 
and for the short-time regime such that the correlation (see equations ([5]) and (|6])) can also be worked out by using 
length £(t) ~ t x / z is smaller than the typical lattice side equation (JSJ) and one gets 
(£ < L,M) but slightly larger than the lattice spacing, 
one has that for k — 1 equation © becomes 

m(t) ~ t\ (9) 

and 

where the exponent 9 — Xa ~^l v is the scaling exponent of 



U 2 {t)~t d '\ (14) 



U±(t) ~ t d l\ (15) 



the initial increase of the magnetization and equation ^ 
holds for t x °/ z < 1[13. 

Furthermore, for e = and mo = the scaling be- respectively, 
havior of the fluctuations of the second moment of the Summing up, by properly measuring the dynamic be- 

magnetization, which give the susceptibility according to havior of the system at criticality one should be able to 

equation Q, is given by[H] determine the exponents, 6 (equation ©) and (d—2/3/u) 



z 



x (t) ~tW- 2 /»/")/*. (10) 

It is worth mentioning that if the hyperscaling relationship 
vd — 2/3 = 7 holds, one has the well-known dependence 
X[t) - Vl vz , where 7 is the susceptibility exponent. 

On the other hand, by starting from a well ordered 
initial configuration, e.g., for mo = 1, the scaling approach 
given by equation ((5J) for k — 1 becomes 

m{e,t)r,t-^ vz m{t^ vz e). (11) 



(equation (jTUJ) ) , by starting from disordered initial config- 
urations; and fi/vz (equation (jlllO . 1/vz (equation (|13[1 ). 
and d/z (equation ((H])), by starting from ordered initial 
configurations. 

Usually the dimensionality of the system is known so 
that this set of five independent determinations becomes 
redundant in order to evaluate four critical exponents. 
However, as will be shown later on, in the case of a system 
in a thermal gradient the complete set is essential for the 



Therefore, just at criticality (e = 0), one should observe a determination of the effective dimensionality involved in 

power-law decrease of the initial magnetization according the scaling relationships. Of course, in our gradient simu- 

to|18j lations we obtained (simultaneously and by means of sin- 

m(t) ~ t~^l vz ', e = 0, (12) gle sets of simulations) information over a wide range of T 
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(Ti < T c < T2), while the above discussed scaling relation- where v is the correlation length exponent, 
ships are expected to hold only close to the critical region, 



i.e., for columns of coordinates such that x r ~ < & — ^4 L 
(see equation JT])). 

On the other hand, equation © is also useful in order 
to describe the stationary scaling regime. In fact, for t — > 
00 the value of mo becomes irrelevant and one gets the 
standard scaling relationship for the magnetization given 
by[H] 

m{e,L) = L-^ l 'fh{eL 1 ' v ), (16) 

where m is a suitable scaling function. Furthermore, just 
at criticality (e = 0) equation (|16p yields 



4 RESULTS AND DISCUSSION 



Let us begin our presentation with results correspond- 
ing to dynamic measurements obtained by starting the 
simulations with fully ordered configurations such that 
m(i,t = 0) = 1, Vi, which corresponds to the ground 
state configuration at T = 0. By following the standard 
procedure the sample is suddenly quenched to the temper- 
ature of the heat bath with the proper thermal gradient. 
So, each column of the sample evolves, at its own teni- 
ae = 0) ~ L , (17) perature, towards its stationary configuration (t —> 00), 
which reflects the fact that the magnetization vanishes at and in particular we focus our attention on the critical re- 
criticality in the thermodynamic limit (L ->■ 00). S ion - According to equation ©, or cquivalently to equa- 
Also in the stationary limit, the spin-spin correlation tion G9, one expects to measure a power-law behavior 
function (equation ©) has its own scaling relationship, of thc Physical observables just at criticality, while devi- 
given bv[T{?] ations from that type of behavior have to be found away 
g(e — r) ~ 7-~( d ~ 2 + T ') (18) from T c . By plotting all the relevant measured observ- 

. . . . . . , ables (m, dl Z^ , U2 and UA as obtained for samples of 

where rj is the critical exponent associated with spatial a(e > 



different sizes (100 < L < 1000), the best set of power 
laws is obtained for T c = 2.2691(3) 20], which is in full 
agreement with the Onsager critical temperature|15j of 
the Ising model in d = 2. Figure 2 shows the time de- 
pendence of the second-order cumulant (see equations ([5]) 
and |H])) as obtained close to criticality by using a square 

lattice of side L = 1000. Here, the thermal gradient is set 

where & is the correlation length, at temperature Ti , which , , ™ „ „ Jrri „ . ,, 

s ' between T± — 2.0 and T 2 — 2.4, while the results . 



correlations, which is related to other exponents by the 
relationship 7 = u(2 — r/)[T5]. 

Furthermore, for homogeneous systems away from the 
critical temperature, the spin-spin correlation function ex- 
hibits an exponential decay behavior given by: 

g(Ti,r)~e- r K<, (19) 



are aver- 



diverges close to criticality, according to , AT innnnj-a- , ^ u A c 4 c 

•" aged over Nc — 10000 different runs. lrom the best ht of 

& ~ A \Ti - T c \~ u , (20) the data shown in figure 2 we obtained d/z = 0.474(6) gD] . 
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Fig. 2. (Color online) Log- log plot of the second-order cumu- 
lant versus time as obtained for a lattice of side L — 1000 and 
averaging over Nc = 10000 different runs. The full line corre- 
sponds to the best power-law behavior that is identified as the 
critical point. Notice that dashed (dashed-dotted) lines show 
upward (downward) deviations that are typical for T > T c 
(T < T c ). The best fit of the solid line yields d/z = 0.474(6). 
More details in the text. 

Our measurements obtained by using initially ordered 
configurations are completed by log-log plots of the time 
dependence of the magnetization and its logarithmic deriva- 
tive, as shown in figures 3(a) and 3(b), respectively. The 
best plots of the lines, obtained for data corresponding 
to the same column that was identified with the critical 
temperature in figure 2, yield P/vz — 0.058(1) 20], and 
l/vz = 0.4857(4) [ID], respectively. 



«1 




Fig. 3. (a) (Upper panel) and (b) (lower panel) show log-log 
plots of the magnetization and its logarithmic derivative ver- 
sus time, respectively. Results obtained at criticality by using 
samples of side L — 1000 and by averaging over Nc = 10000 
different runs. In both panels the straight lines correspond to 
the best fits of the data shown by means of circles. The lin- 
ear gradient of the thermal bath is set between T\ = 2.0 and 
T 2 = 2.4. More details in the text. 

librium conditions is identified with the magnetic suscep- 
tibility, as shown in figure 4. 



Furthermore, initially disordered configurations with 
On the other hand, by starting the simulations with a preset (vanishing) initial magnetization are suitable for 
fully disordered initial configurations such that mo = 0, the measurement of the initial increase of the order param- 
which corresponds to samples in contact with a thermal eter, which after the scaling arguments discussed within 
bath at T — oo that are then suddenly annealed to a de- the context of equation ©, is expected to be observed 
sired thermal gradient, we measured the time evolution of for times such that t << m z ^ x ° . Figure 5 shows log-log 



the fluctuations of the order parameter that under equi- plots of m versus t as obtained for the different values of 
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Fig. 4. Log-log plot of the fluctuations of the order parameter 
versus t, as obtained for the column identified with the critical 
temperature in the plots shown in previous figures. Data corre- 
sponding to L = 1000 and taking a linear gradient set between 
Ti = 2.0 and T 2 = 2.4. Averages are taken over N c = 10000 
different initial configurations. The best fit of the data, shown 
as a full line, has slope y/vz — 0.3570(4). 

mo (0.035 < tuq < 0.050). The data shown in figure 5 
correspond to the previously identified "critical" column. 
The best fits of the curves give estimates of the initial 
increase exponent 9 that, after a proper extrapolation to 
the limit tuq (not shown here for the sake of space), 
yields 6 = 0.196(6)[D]. 

Based on the results obtained by means of dynamic 
measurements only, we are in a condition to outline a few 
interesting preliminary conclusions on the critical behav- 
ior of the Ising magnet in contact with a gradient thermal 
bath. On the one hand, the best fits of all physical observ- 
ables were found at the same column (i.e., the same tem- 
perature) of the sample, which has been identified as the 
critical temperature T c = 2.2691 (3) 20J . This value is in 
full agreement with the well known exact result early eval- 




Fig. 5. Log-log plot of the magnetization versus t, as ob- 
tained for the column identified with the critical temperature 
in the plots shown in previous figures. Data corresponding to 
L — 1000 and taking a linear gradient set between Ti = 2.0 and 
T-2 — 2.4. Averages are taken over Nc ~ 10000 different ini- 
tial configurations. Each curve corresponds to a different initial 
magnetization mo, as indicated, and the dashed lines are the 
fits of the data. The full straight line corresponds to 8 = 0.196, 
i.e. the value of the initial increase exponent obtained by ex- 
trapolation to mo — ¥ 0. More details in the text. 

uated by Onsager, i.e., Tonsager — 2.2692. On the other 
hand, by using fi/vz = 0.058(1) and l/vz = 0.4857(4)[20] 
as determined by means of the measurements of the mag- 
netization and its logarithmic derivative, respectively, we 
obtained (3 — 0.120(2)[5D], in excellent agreement with the 
exact value f3 = 1/8 = 0.1 25 [15]. Additionally, by assum- 
ing that the hyperscaling relationship dv — 2/3 — 7 holds, 
just dividing by vz one can write 

^-^=0, (21) 

z vz vz 

where equation (|2ip can be interpreted as a "dynamical" 
hyperscaling relationship. Then, by replacing the mea- 
sured exponents in equation (I2T1) wc obtain 0.001(9) on 
the right-hand of the equality, a result that strongly sup- 
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ports the validity of hyperscaling. On the other hand, we 
found that the exponent for the initial increase of the order 
parameter (0 — 0.196(6)) is in full agreement with previ- 
ous numerical results obtained by applying the Metropolis 
dynamics to the two-dimensional Ising magnet (d = 2) in 
a homogeneous bath, namely, 9 = 0.191(1)[2T]. However, 
it should be mentioned that this exponent depends on the 
dynamics used (e.g. Metropolis, Glauber, Heat-Bath, etc.) 
and that our Gradient-Metropolis dynamics may not give 
the same exponents as the proper standard Metropolis 
dynamics. 

Let us now point our attention to stationary results 
obtained after disregarding Nd = 5 x 10 5 Monte Carlo 
steps, in order to allow for the stabilization of the sample, 
and evaluated during the subsequent time interval of 5 x 
10 5 Monte Carlo steps. 

Figure 6(a) shows plots of magnetization profiles (i.e., 
plots of 77i (i) versus Ti, where Tj is the temperature of the 
ith column), as obtained for samples of different sizes. 
Since the temperatures at the extremes of the sample 
are kept constant at T\ = 0.64 and T2 = 3.40, each 
curve in figure 6(a) corresponds to different gradients. It 
follows that finite-size effects are negligible in the low- 
temperature range, e.g., for T < 2.1, while these effects 
become slightly evident above criticality. 

It is worth mentioning that the results shown in fig- 
ure 6(a) can be replotted in order to test whether the 
standard scaling also holds for gradient measurements. In 
fact, figure 6(b) shows log-log plots of m(Tj, L)L^/ V versus 
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Fig. 6. (Color online) (a) Plot of the magnetization profiles 
versus T as obtained for samples of different thermal gradients. 
The temperatures at the extremes of the sample are T\ — 0.64 
and T2 = 3.40, and results are averaged over Nd = 5 x 10 5 
Monte Carlo steps after disregarding the initial 5 x 10 5 Monte 
Carlo steps in order to allow for the stabilization of the sam- 
ples, (b) Scaling plots of the data already shown in (a) as ob- 
tained by using equation (|16p . More details in the text. 



\Ti -T C \L X I V as obtained by using the data already shown 
in figure 6(a). The quality of the collapse supports the va- 
lidity of the standard scaling Ansatz for the case of our 
gradient measurements; however, the terms of correction 
to scaling, which are neglected in our analysis, cannot be 
disregarded. On the other hand, just by selecting the crit- 
ical column one has that equation (|17[) simply gives the 
monotonic decay of the magnetization when the system 
size is increased, namely m(T — T c ) ~ L~PI V , a result 
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Fig. 7. Log-log plot of the magnetization, measured close to 
criticality, versus the system side L. Data corresponding to the 
critical column (open circles), and two adjacent columns: T < 
T c , open squares, and T > T c , open diamonds, respectively. 
Results obtained by averaging over 5 x 10 5 Monte Carlo steps 
after disregarding the first Nn = 5 x 10 5 Monte Carlo steps 
in order to allow for the stabilization of the system. The full 
line has slope —/3/v — 0.125(9), and corresponds to the best 
fit of the data. The other two lines correspond to the fits of 
the data from adjacent columns and give the error bars to the 
result 20 . More details in the text. 
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Fig. 8. (Color online) Plots of the second- (a) and fourth- 
order (b) cumulants U2 and Ua versus the temperature of the 
column, as obtained for different values of the system side L as 
indicated. The vertical straight lines on each graph indicate the 
intersection points and are shown for the sake of comparison. 



that has been verified in figure 7 where the best fit of the 
data corresponds to fi/v = 0.125(9)[20]. 

As already established in the field of critical phenom- 
ena, the so-called cumulants (see equations (O and ©) 
are suitable functions of the moments of the order param- 
eter distribution function whose pre-scaling factor, disre- 
garding high-order finite-size scaling corrections, is inde- 
pendent of the system size. So, plots of the cumulants 
versus the control parameter (i.e., the temperature of the 
gradient thermal bath) should exhibit a common inter- 
section point, as it is indicated with full straight vertical 



lines in the data shown in figures 8(a) and 8(b), obtained 
with our gradient system for {/ 2 (L, T) and U&(L,T), re- 
spectively. A careful inspection of the critical region close 
to the interaction points reveals very small but system- 
atic shifts of the intersection points between curves cor- 
responding to adjacent system sizes. After proper extrap- 
olation of the intersection points to the thermodynamic 
limit (not shown here for the sake of space), we obtain 
T c (oo) = 2.284(1) [20, and T c (oo) = 2.283(3) [20], for data 
corresponding to Ui and C/4, respectively. 
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On the other hand, plots of the fluctuations of the or- 
der parameter, which are identified with the susceptibility 
in standard measurements, show peaks close to criticality 
as expected (see figure 9). In fact, it is known that the 
susceptibility exhibits rounding and shifting effects as the 
size-dependent "critical" temperature of the peaks (T C (L)) 
converges toward the true critical temperature according 
to a standard finite-size scaling relationship given bv|18] 

T C (L) = T c (oo) + constant L~ 1/v '. (22) 

So, the inset of figure 9 shows plots of T C (L) versus L _1 
(here we assume v — 1 for the correlation length exponent 
as justified below) that yield T c (oo) = 2.296(1) [213 - 

In additional simulations, we studied the spin-spin cor- 
relation functions, both parallel and perpendicular to the 
gradient direction. We measured the correlation function 
in the perpendicular direction, by using equation (0 for 
different temperatures, both lower and higher than T c ~ 
2.26918. We used a rectangular lattice of size L = 300, 
M = 4000, and set the temperatures of the extremes of 
the sample at T\ = 0.8 and Ti = 3.74, respectively. For 
each measurement corresponding to a certain temperature 
outside criticality, we plotted the correlation function ver- 
sus the distance r, in log-linear plots as shown in figure 
10(a). Then, by fitting the data to an exponential decay, 
according to equation (IT91 , we obtained the values of the 
correlation length £j(|Tj — T c \). The inset of figure 10(a) 
shows a linear-linear plot of £(T) versus 1/(T — T c ) that 
yields a straight line. This behavior confirms that equa- 
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Fig. 9. (Color online) Plots of the fluctuations of the mag- 
netization versus the temperature of the column, obtained by 
means of simulations with systems of different side, as listed in 
the figure. The inset shows the extrapolation of the effective 
"critical" temperature according to equation (|22[) . which yields 
T c (oc) = 2.296(1) as the best fit (open circles). Notice that ex- 
trapolation of data corresponding to the columns adjacent to 
the critical one are also included, such that open diamonds and 
open squares correspond to T < T c and T > T c , respectively. 
More details in the text. 

tion ([2T)f holds with v — 1, and from the slope we also 
obtained the prefactor given by A = 1.16(4). 

On the other hand, for the correlation function corre- 
sponding to the column of the system with temperature 
closest to T c ~ 2.26918, we performed a power-law fit ver- 
sus the distance r, according to equation (JT5J) (see figure 
10(b)), which yields d-2 + n = 0.238(3). By eliminating 
the value of d from the exponent equation with the use 
of relationship dv — 2/3 — 7, and by using values already 
calculated from both stationary and dynamic simulations, 
we obtained n — 1.22(1). This value is in excellent agree- 
ment with n — 5/4 = 1.25, found by taking d = 1 in 
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Fig. 10. (Color online) (a) Log-linear plots of the vertical spin- 
spin correlation function g(T, r) versus the distance r. Different 
curves correspond to different temperatures away from critical- 
ity, as indicated. For each curve an exponential decay fit was 
made, obtaining values of the correlation length f(T). The in- 
set in (a) shows the values of £ versus 1/\T — T c \. The linear 
fit indicates that equation (|20p holds, with an exponent v = 1. 
Furthermore, the slope of the fit gives us the value of the pre 
factor A in equation ([20]), which is found to be A = 1.16(4). 
(b) Log-log plot of g(T,r) versus r, as obtained for the critical 
temperature, with its corresponding power-law fit (full line). 
The exponent obtained corresponds to d — 2 + r\ = 0.238(3). A 
rectangular lattice was used, with L — 300 and M — 4000. The 
temperature gradient was set between Ti — 0.8 and T2 = 3.74. 
Results obtained by averaging over 5 x 10 s Monte Carlo steps 
after disregarding the first Nd = 5 x 10 s Monte Carlo steps in 
order to allow for the stabilization of the system. More details 
in the text. 
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the relationship that results from combining 7 = v(2 — rf) 
and dv — 2/3 = 7 with the exact exponent values of the 
two-dimensional Ising Model, i.e., v = 1 and j3 = 1/8 [TSj. 

On the other hand, the correlation function in the di- 
rection parallel to the thermal gradient has also been mea- 
sured, by using a system of L = 5000 and M = 100, with 
a gradient between temperatures T± — 0.8 and Ti = 3.74, 
respectively. The correlation function studied corresponds 
to g(T,r) at the critical temperature, namely, g(T c ,r). It 
is worth mentioning that correlations along the direction 
of the applied gradient involve spins at different tempera- 
tures, so a quantitative study such as the one done for ver- 
tical correlations is no longer possible. Instead, we gain an 
insight into the qualitative behavior of the Ising ferromag- 
nct in a thermal gradient. In fact, figure 11 shows g(T c , r) 
versus the distance r. For r > 0, we found that correla- 
tions between the spins in the critical column and those 
placed in columns at higher temperatures are almost neg- 
ligible. On the other hand, for r < correlations between 
the spins located in the critical column and those placed 
in columns at lower temperatures are measured. Here, the 
correlations first abruptly decrease and then they become 
almost constant for colder temperatures. The inset in fig- 
ure 11 shows a log- log plot of a zoomed view of the area 
near criticality, where absolute values of r were taken, in 
order to account for both branches of g(r), namely, for 
r > and r < 0, respectively. The full line is a plot 
of a power-law decay with exponent 1/4 = 0.25 and is 
shown for the sake of comparison. By comparing the mea- 
sured points with the full line, it can be inferred that close 
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Fig. 11. Plots of the correlation function g(T,r), where the 
origin is taken just at the critical temperature. Here, for r < 
(r > 0) one has T < T c (T > T c ). A rectangular lattice was 
used, with L — 5000 and M = 100. The temperature gradient 
was set between Ti = 0.8 and T2 = 3.74, respectively. Results 
obtained by averaging over 5 x 10 5 Monte Carlo steps after 
disregarding the first Nd = 5 x 10 5 Monte Carlo steps in order 
to allow for the stabilization of the system. The inset shows an 
enlarged view of the area close to criticality, on log-log scale. 
Squares (circles) correspond to correlations with spins in the 
direction of T < T c (T > T c ). The full line is a plot of the 
power law of r with exponent 1/4 = 0.25 and is shown for the 
sake of comparison. Absolute values of r have been taken, in 
order to show both branches of g(r). More details in the text. 

enough to the critical column correlations have a func- 
tional behavior in agreement with the well known scaling 
law r~( d ~ 2+ ''). This result reflects the fact that for the 
considered gradient ^ = 0.588 x 10~ 3 , the region close 
to r = is almost at the critical temperature. 

Let us now discuss the results obtained in the present 
paper and perform comparisons to existing numerical re- 
sults and/or exact values corresponding to the Ising mag- 



Type 


Observable 


Tc 


Dynamic 


Magnetization 


2.2691(3) 


Dynamic 


Cumulant U2 


2.2691(3) 


Dynamic 


2 nd order magnetization 


2.2691(3) 


Stationary 


Cumulant U2 


2.284(1) 


Stationary 


Cumulant U4 


2.283(3) 


Stationary 


Scaling of m 


2.282(3) 


Stationary 


Extrapolation of T C (L) from x 


2.296(1) 


Exact [15] 




~ 2.26918 



Table 1. List of critical temperatures (third column) obtained 
by means of dynamic and stationary measurements, as listed 
in the first column. More details in the text. 



net in a homogeneous thermal bath. Table I summarizes 
over seven independent estimations of the critical temper- 
ature as follows from the measurement of different phys- 
ical observables in our gradient Ising system. In all cases 
we obtained an excellent agreement with the exact value, 
pointing out that the proposed gradient method is suitable 
for accurate determinations of critical points. 

Pointing now our attention to the measured critical 
exponents, we have already shown that our dynamic mea- 
surements of d/z, /3/vz, and 7 jvz are fully consistent with 
a sort of "dynamic" hyperscaling relationship (equation 
(|21jl ). However, since the dimensionality entering in equa- 
tion (|21[) appears in a quotient between exponents, more 
precisely as d/z as follows from dynamic measurements 
of the cumulant, one needs independent measurements in 
order to obtain an unbiased estimation of d. The same rea- 



16 




n it i • l A n mi t ■ 

Mugha and Albano: lhe ism 


Exponent 


Present work 


Exact/numerical value 


P 


0.120(2) 


1/8=0.125 


V 


0.96(7) 


1 


V 


1.238(3) 


5/4=1.25 (Details in the text) 


7 


0.736(1) 


3/4=0.75 (Details in the text) 


z 


2.16(4) 


2. 166(71122] 


e 


0.196(6) 


0.191(1)^11 


d 


1.02(8) 





Table 2. List of the critical exponents (second column), along 
with the value of the dimensionality, obtained by combining 
dynamic and stationary measurements. The third column lists 
the exact value, in the first four rows, and the comparison value 
from numerical simulations, in the fifth and sixth rows. 



soning applies to the dynamic exponent z. Therefore, we 
conclude that in order to obtain all the critical exponents 
we have to combine measurements corresponding to both 
the short-time dynamic simulations, namely, d/z, (3/vz, 
and 7/^2:, and stationary simulations, namely, f3/v and 
d — 2 + rj. The results are shown in Table II, along with 
the values for the exponents that are known exactly (/3, v, 
rj, and 7), and the best available numerical results (z and 
9), for the two-dimensional Ising model with a homoge- 
neous thermal bath. It is worth mentioning that the error 
bars in each value were propagated from those of the orig- 
inal exponents, which in turn were calculated considering 
the values of the physical observables from the adjacent 
columns (at different temperatures) to the critical one [20]. 
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There is an interesting result that involves both 7 and 
77 because the values obtained are different in a quantity of 
1 as compared with the exact values for the Ising ferromag- 
net for d — 2, namely 7 = 7/4 and 77 = 1/4. Furthermore, 
if we determine the values of 7 and r\ from the scaling rela- 
tionships dis — 2/3 — 7 and 7 = i^(2 — 77) by using the values 
of the exponents found in this work, we roughly reach the 
same results, namely 7 = 0.74(8) and -q = 1.23(8) re- 
spectively. The reason for this apparent disagreement be- 
comes evident when we calculate, from the combination 
of dynamic and stationary results, the value of the effec- 
tive dimensionality, obtaining d = 1.02(8), which strongly 
suggests d = 1 as judged by the uncertainty interval. 
If we now use d = 1 and re-calculate the exponents of 
the Ising ferromagnet, we obtain 7 = 3/4 = 0.75 and 
rj = 5/4 = 1.25, which are in perfect agreement with 
the exponents found by means of our simulations. This 
value for the effective dimensionality d takes into account 
the fact that due to the presence of a thermal gradient, 
the thermodynamic functions are calculated as profiles of 
columns, thus lowering the effective dimensionality exactly 
by one unit. 

5 CONCLUSIONS 

We have studied, by means of Monte Carlo numerical sim- 
ulations, the dynamic and stationary critical behavior of 
the two-dimensional Ising ferromagnet, in contact with a 
thermal bath that exhibits a linear gradient between low- 
and high- temperature extremes at T± and Ti {T2 > T\), 
respectively. It is worth mentioning that our proposed ap- 
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proach, is based on the assumption that the ferromagnet 
is in thermal equilibrium with the heat bath that intrinsi- 
cally has a thermal gradient. This is in contrast to previous 
studies by other authors [l][2l[3j|4] who placed only the ex- 
tremes of the magnet in contact with two thermal baths at 
T\ and T2, and subsequently measured the thermal con- 
duction between these baths. While for the numerical im- 
plementation of this "conductivity approach" one needs 
to use a suitable "Creutz demon" algorithm [TO], our more 
straightforward approach can be implemented by means of 
standard algorithms, e.g., the Metropolis dynamics, just 
as in the present paper. Of course, our whole sample is out 
of equilibrium, but a key feature is that each column of the 
magnet (or sample in general) can be considered in equi- 
librium with the gradient thermal bath that is in contact 
with it. In this way, now one not only has information of 
the system under study for a wide range of temperatures 
(Ti < T < T2), but also by means of a suitable choice 
of the temperatures of the extremes of the sample, such 
that T\ < T c < T2, where T c is the critical temperature, 
one can deal with samples simultaneously exhibiting the 
ordered and the disordered phases. In this way, instead 
of measuring average values of the physical observables at 
each T, one actually measures "thermal profiles" of the 
observables averaged over (d — l)-dimensional columns of 
the sample. For this reason, the effective dimensionality 
entering in the scaling relationships is d = 1. 

Focusing our attention on the dynamic measurements, 
the measured thermal profiles allow us to simultaneously 
follow the time evolution of the sample for all tempera- 
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tures of interest and, by drawing suitable plots, quickly 
identify the critical temperature (corresponding to a par- 
ticular column) and determine the relevant critical expo- 
nents. Similarly, for the case of stationary measurements, 
we take advantage of the wide information available (criti- 
cal temperature, critical exponents, correlation functions, 
etc.) that can be obtained just by performing a single sim- 
ulation. The results obtained for the critical temperature 
and critical exponents are summarized and discussed in 
the context of Tables I and II, respectively. Based on that 
analysis, we conclude that the figures are remarkably accu- 
rate when one considers the computational effort involved. 
In fact, full agreement is obtained when our results are 
compared with exactly known values, e.g., for the critical 
temperature T c , and the critical exponents /?, v, 7, and 77, 
as well as with the best available data taken from numer- 
ical simulations of other authors, e.g., for the exponents 9 
and z. 



Based on the results presented and discussed in this 
paper, we conclude that the study of critical systems in 
the presence of a gradient of the suitable control parameter 
is a useful and powerful tool to gain rapid inshight on the 
behavior of the system, which could eventually be studied 
in further details by using more sophisticated methods . 
Furthermore, the studies of material systems under gra- 
dient conditions could shed light on the understanding of 
a wide variety of experimental and technologically useful 
situations. 
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